
#########################
#                       #
#     Simulations       #
#                       #
#########################
rm()
memory.limit(20000)
library(plm)
library(knitr)
library(survey)
library(matrixStats)
library(nlWaldTest)

#Create table for results
table.ARDL<- data.frame(`Alpha1_SRE=1`="alpha SRE = 1", `ARDL_Beta1_Type2`="Beta1 Type 2", `ARDL_LRE_Type1_SRE=1`="LRE Type 1 SRE = 1")

rownames(table.ARDL)[1] <- "T"

for (k in 5:100) {
  #First, define the parameters used by the data-generating process:
  #T is the number of waves of data including the first wave needed only to produce the lag dependent variable for the second wave. 
  
  rho <- 0.8
  sig1 <- 1 
  sig2 <- 1 
  alpha1 = 1
  beta1 = 1
  beta2 = -1
  LT <- k*5
  T<-LT+50
  reps <- 100000
  
  set.seed(123)
  
  #The following function generates a synthetic dataset of desired dimensions (T time points) and distribution parameters:
  
  generate <- function(T, rho, sig1, sig2) {
    y <- matrix(0, T)
    xI <- matrix(0, T)
    xS <- matrix(0, T)
    for (t in 1:T) {
      xxS <- if (t>1) xS[t-1] else 0
      xS[t] <- rho * xxS + rnorm(1, sd = sig1)
      xxI <- if (t>1) xI[t-1] else 0
      xI[t] <- xxI + rnorm(1, sd = sig2)
      y[t]<- xI[t] + xS[t]
    }
    data.frame(t = seq(1,T-50),
               x1=c(xS[(51):(T)]),
               x1_lag=c(xS[(50):(T-1)]),
               y = c(y[(51):(T)]),
               y_lag = c(y[(50):(T-1)]))
  }
  
  #Now we generate datasets T time points 
  
  ds <- replicate(n = reps,
                  generate(T = T,
                           rho = rho,
                           sig1 = sig1,
                           sig2 = sig2),
                  simplify = FALSE)
  
  
  #########################
  #                       #
  #   OLS-ARDL            #
  #                       #
  #########################
  
  #Now we fit the ARDL model to the datasets and save the results
  set.seed(421)
  lms<-lapply(ds,
              function(d) {
                lm(y~y_lag + x1 + x1_lag,
                     data = d)
              })
  
  #Let's check the sampled lm parameters:
  est_param <- sapply(lms, coef)
  est_confint <- sapply(lms, function(lm.model){confint(lm.model, parm=c(2,3,4))})
  lr_confint<-lapply(lms, function(lm_model){nlConfint(lm_model, "(a[3] + a[4])/(1-a[2])")})
  power_alpha1<-ifelse(est_confint[4,]<0 | est_confint[1,]>0 , 1, 0)
  power_beta1<-ifelse(est_confint[5,]<0 | est_confint[2,]>0 , 1, 0)
  power_beta2<-ifelse(est_confint[6,]<0 | est_confint[3,]>0 , 1, 0)
  power_lrun<-sapply(lr_confint, function(lr_confint) {ifelse(lr_confint[2]>0 | lr_confint[3]<0, 1, 0)})
  power<-c(mean(power_alpha1), mean(power_beta1), mean(power_beta2), mean(power_lrun))
  
  table.ARDL[k-3,]<- data.frame(`Alpha1_SRE=1` = mean(est_param[2,]),
                                `ARDL_Beta1_Type2` = 1 - power[2],
                                `ARDL_LRE_Type1_SRE=1` = power[4])
  
  rownames(table.ARDL)[k-3] <- (k*5)


}
write.csv(table.ARDL, "ARDL--SRE1.csv")

